Notas previas:
Los objetivos específicos de la práctica se dividen en varias fases analíticas:
El análisis metagenómico explora cómo la microbiota inicial puede influir en la respuesta del sistema nervioso central a la inflamación sistémica, buscando correlaciones entre la microbiota y las respuestas a la sepsis.
El análisis de RNAseq busca identificar mecanismos moleculares y celulares que contribuyen a la neurodegeneración y neuroinflamación en un modelo murino de sepsis, con un enfoque en las diferencias de género.
El análisis de scRNAseq profundizará en la respuesta celular específica a la sepsis, centrando la atención en células como los macrófagos microgliales y células linfoides para entender las dinámicas de población celular.
Finalmente, en futuros estudios se planea una integración multiómica de los datos de RNAseq, scRNAseq, y metagenómica para descubrir interacciones complejas y posibles dianas terapéuticas, utilizando datos proporcionados por el grupo de Oftalmología Experimental de la Universidad de Murcia.
Se va a realizar un análisis comparativo del metagenoma de ratones criados en la Universidad de Murcia frente a aquellos criados en la Universidad de Leuven, Bélgica. Este estudio se motiva por observaciones preliminares que indican que los ratones de Murcia exhiben una mayor sensibilidad y una respuesta inflamatoria exacerbada al tratamiento con LPS (lipopolisacárido), a diferencia de los ratones de Leuven, que muestran una reacción atenuada o no lo presentan.
La secuenciación metagenómica de estas muestras se llevó a cabo utilizando la tecnología Ion Torrent. Para el análisis y visualización de los datos obtenidos, se empleará QIIME 2 (Quantitative Insights Into Microbial Ecology 2), una herramienta de software de código abierto especialmente diseñada para estudiar la ecología microbiana a partir de secuencias de ADN. Este enfoque permite identificar diferencias en la composición de la microbiota que podrían explicar las variaciones en la respuesta al LPS entre las cohortes murinas de ambas universidades.
El primer paso crucial en el uso de QIIME 2 es la importación
adecuada de los archivos FASTQ. Para nuestro estudio, disponemos de 19
archivos FASTQ correspondientes a 19 muestras diferentes, que ya están
demultiplexadas, según se detalla en el archivo
metadata.tsv. Para la importación, utilizamos el formato
Casava 1.8, y optamos por la estrategia de importación mediante un
archivo de manifiesto. Este archivo de texto simple es fundamental para
que QIIME 2 pueda interpretar correctamente la estructura de los datos
de entrada. En él se mapean las muestras a sus respectivos archivos de
secuencias, especificando la ubicación exacta de los archivos y la
manera en que deben ser procesados. El manifiesto debe incluir, al
menos, dos columnas: una para el identificador de cada muestra y otra
para las rutas a los archivos de secuencias (consultar
manifest.csv).
Una consideración esencial al importar archivos Casava 1.8 mediante manifiesto en QIIME2 es que los nombres de los archivos no deben contener el carácter “_”. Este carácter puede causar problemas durante el procesamiento posterior, ya que QIIME 2 podría no reconocer el archivo correctamente.
El manifiesto se genera a partir de las siguiente linea de comandos:
echo "sample-id,absolute-filepath,direction" > manifest.csv
ls *.fastq | awk -v pwd="$PWD/" '{sub(/\.fastq$/, "", $1); print $1 "," pwd $1 ".fastq,forward"}' >> manifest.csvLo siguiente es ejecutar la línea de comandos para importar los ficheros FASTQ a qiime para el procesamiento de los mismos. Para visualizar la calidad de las lecturas de los ficheros FASTQ se ejecuta seguidamente la siguiente linea:
qiime tools import \
--type 'SampleData[SequencesWithQuality]' \
--input-path manifest.csv \
--output-path demux.qza \
--input-format SingleEndFastqManifestPhred33
qiime demux summarize \
--i-data demux.qza \
--o-visualization demux.qzvLa Figura 1 de la izquierda muestra la calidad promedio de las lecturas de nuestras muestras. Observamos que la calidad se mantiene en un promedio de 30 para las primeras 280 bases, pero luego se vuelve inestable. Por esta razón, decidimos excluir esta sección inestable de las secuencias.
Dado que nuestra herramienta está configurada para procesar lecturas de secuenciación de Illumina, necesitamos realizar algunos ajustes en el protocolo de procesamiento. Conforme a Catozzi et al. 2019, es recomendable recortar las primeras 15-20 bases de cada secuencia para eliminar los primers utilizados en la secuenciación. Además, seguimos la recomendación de Salipante et al., 2014, de descartar lecturas con menos de 100 bases.
Para realizar estos ajustes, empleamos un workflow en Galaxy. Iniciamos con la herramienta “FASTQ Trimmer by column” (Galaxy Version 1.1.5) para eliminar las 20 primeras bases y recortar las secuencias hasta la base 280. Posteriormente, utilizamos “Filter FASTQ reads by quality score and length” (Galaxy Version 1.1.5) para filtrar aquellas secuencias que superen las 100 bases de longitud. Finalmente, estas secuencias filtradas se importan de nuevo a QIIME2 para su análisis posterior, como se ilustra en la Figura 1B.
Existen dos estrategias para calcular la diversidad de una muestra. Mediante Unidades Taxonómicas Operativas (OTUs) o mediante inferencia de Variantes de Secuencia de Amplicones (ASVs). Las OTUs se agrupan por similitud de secuencia, generalmente con un umbral del 97%, lo que significa que las secuencias que son al menos un 97% similares se consideran la misma OTU. Las ASVs no se agrupan por similitud; cada variante única de secuencia se trata como una entidad distinta.
Según el manual de deblur de QIIME2, que genera OTUs, este comando está diseñado específicamente para datos generados a partir de protocolos de amplicones 16S en plataformas Illumina: “This mode of operation should only be used when data were generated from a 16S amplicon protocol on an Illumina platform”. Siguiendo las recomendaciones de Catozzi et al. y diversas fuentes en blogs, se aconseja emplear la herramienta DADA2 para este proceso. A diferencia del método tradicional que genera OTUs, DADA2 permite la creación de ASVs, que proporcionan una resolución taxonómica más precisa. Para implementar esto, ejecutamos lo siguiente:
qiime dada2 denoise-single \
--i-demultiplexed-seqs demux.qza \
--p-trunc-len 0 \
--p-n-threads 0 \
--o-representative-sequences rep-seqs-dada.qza \
--o-table table-dada.qza \
--o-denoising-stats denoising-stats.qza| Statistic | count | min | max | mean | range | std |
|---|---|---|---|---|---|---|
| Value | 6117 | 25 | 260 | 230.386 | 235 | 34.2403 |
Vemos que se han inferido un total de 6117 secuencias distintas con un tamaño medio de 230 pb (Tabla 1).
Para elegir un sampling-depth adecuado, primero necesitamos examinar la distribución de las profundidades de secuenciación de nuestras muestras. Generalmente hay que elegir una profundidad que mantenga el máximo número de muestras posibles mientras se sacrifica lo menos posible en términos de número de secuencias. La mediana de la frecuencia por muestra que encontramos que es de 101.820 features. Sin embargo, decidimos escoger una sampling depth de 30.000 features, recogiendo así el 25,25% de las features en el 78,95% de las muestras.
La generación de un arbol filogenético es fundamental para entender las relaciones evolutivas entre las secuencias de especies en nuestras muestras. La generación de un árbol filogenético se basa en diferentes métricas como la diversidad filogenética de Faith, UniFrac ponderado y no ponderado. El método se basa en un alineamiento múltiple de las secuencias para generar la filogenia en base a la similitud de las mismas. Este paso es esencial para el posterior análisis de métricas de diversidad alfa y beta.
Ya tenemos todo lo necesario para hacer el árbol filogenético. Esto se genera mediante los siguientes comandos:
qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree-dada.qza \
--i-table table-dada.qza \
--p-sampling-depth 30000 \
--m-metadata-file metadata.tsv \
--output-dir core-metrics-results \
--p-n-jobs-or-threads auto \
--parallelPara observar si hemos recogido suficiente diversidad en nuestras muestras realizamos un análisis de Alpha rarefacción (Figura 2). Para ello ejecutamos lo siguiente:
qiime diversity alpha-rarefaction \
--i-table table-dada.qza \
--i-phylogeny rooted-tree-dada.qza \
--p-max-depth 30000 \
--m-metadata-file metadata.tsv \
--o-visualization alpha-rarefaction-dada.qzvFigura 2: Gráfico de rarefación. Índice shannon de cada muestra según la profundidad de secuenciación.
El análisis de la gráfica de rarefacción alfa (Figura 2) nos indica que todas las muestras alcanzan un máximo en el índice de Shannon, lo cual sugiere que incrementar el número de secuencias probablemente no revelará un número significativo de especies nuevas. Esto indica que por más profundidad de secuenciación no encontraremos especies nuevas. Además, se observan variaciones en la diversidad alfa entre las diferentes muestras; algunas registran valores más elevados en el índice de Shannon que otras. Estas diferencias podrían atribuirse tanto a variaciones reales en la biodiversidad de las comunidades muestreadas como a la eficacia del muestreo o de los métodos de secuenciación empleados.
En el análisis de metagenómica, es habitual calcular diversos índices de diversidad que proporcionan una medida de la abundancia y riqueza de especies en cada muestra. Estos índices facilitan la comparación de estas características entre diferentes muestras o grupos de muestras.
La diversidad alfa se refiere a la riqueza de especies dentro de una unidad de muestreo específica, que puede ser una muestra individual o un grupo de muestras bajo la misma condición. Para medirla, se utilizan varias métricas como el índice de Shannon, la riqueza de especies observadas, equidad de Simpson, diversidad inversa de Simpson y la diversidad filogenética de Faith, entre otras, que consideran tanto la presencia y ausencia de especies como su abundancia relativa.
Por otro lado, la diversidad beta evalúa la variación o cambio de diversidad entre diferentes ambientes o condiciones, utilizando métricas de distancia que reflejan la disimilitud entre muestras, siendo las más comunes la distancia de Bray-Curtis, el Índice de Jaccard y Unifrac, para determinar cuán similares o diferentes son las muestras entre sí.
Para obtener todas las métricas para estudiar en profundidad ejecutamos la siguiente línea de comandos:
qiime diversity core-metrics-phylogenetic \
--i-phylogeny rooted-tree-dada.qza \
--i-table table-dada.qza \
--p-sampling-depth 30000 \
--m-metadata-file metadata.tsv \
--output-dir core-metrics-results \
--p-n-jobs-or-threads auto \
--parallelCon el código anterior se generan los ficheros para observar estas métricas y poder hacer estudios estadísticos al respecto. Como vemos en la Figura 3 y la Figura 4, en la comparación entre grupos de sexo o entre universidades, respectivamente, no observamos ninguna diferencia en la diversidad o abundancia recogidas en las muestras. Si nos fijamos en profundidad, cada universidad tiene unos outliers (2 en la KU y 3 en la UMU) con una diversidad menor que el resto. Esto podría ser debido a un fallo en la secuenciación o que estas muestras tienen algo distinto al resto de muestras.
Figura 3: Medidas de estudio de diversidad alfa en las muestras distinguiendo por Sexo. De izquierda a derecha: Indice Shannon, Abundancia, Faith e indice de equidad.
Figura 4: Medidas de estudio de diversidad alfa en las muestras distinguiendo por Universidad. De izquierda a derecha: Indice Shannon, Abundancia, Faith e indice de equidad.
Vamos ahora a estudiar la diversidad beta para ver si se parecen las especies de microorganismos entre grupos. Tanto la distancia de Bray-Curtis (Figura 5) como la distancia de Jaccard (Figura 6) son medidas para medir la disimilitud entre muestras. Podemos observar que en el PCoA de ambas distancias se diferencian claramente 4 poblaciones dentro de nuestras muestras y estan claramente definidas por las Universidades, ya que el sexo no es un factor diferencial para separar las poblaciones .
Figura 5: PCA plot de la distancia Bray-Curtis para estipular la beta diversidad entre las muestras. Las muestras de color morado pertenecen a la Universidad de Leuven y las de amarillo a la UMU. Los puntos circulares son muestras de sexo femenino y los cuadrados, masculinos.
Figura 6: PCA plot de la distancia Jaccard para estipular la beta diversidad entre las muestras. Las muestras de color morado pertenecen a la Universidad de Leuven y las de amarillo a la UMU. Los puntos circulares son muestras de sexo femenino y los cuadrados, masculinos.
Si hacemos un estudio estadístico de los datos, observamos que con ambas métricas (Bray-Curtis y Jaccard) hay una diferencia significativa entre las muestras de cada Universidad (Tabla 2 y Tabla 3).
| Group.1 | Group.2 | Sample.size | Permutations | pseudo.F | p.value | q.value | Metric |
|---|---|---|---|---|---|---|---|
| KU | UMU | 15 | 999 | 3.565576 | 0.001 | 0.001 | Bray-Curtis |
| Group.1 | Group.2 | Sample.size | Permutations | pseudo.F | p.value | q.value | Metric |
|---|---|---|---|---|---|---|---|
| KU | UMU | 15 | 999 | 2.801261 | 0.001 | 0.001 | Jaccard |
Cada ASV identificada previamente necesita ser asignada a un taxón específico, y para esta tarea, en QIIME2, se puede emplear un clasificador Bayesiano que ha sido preentrenado. Dicho clasificador se ha derivado de la base de datos Greengenes de OTUs al 99%, la cual ha sido recortada para ajustarse a lecturas de 250 bases correspondientes a la región V4 del ARNr 16S.
Para ello ejecutamos el siguiente código:
qiime feature-classifier classify-sklearn \
--i-classifier gg-2022-10-backbone.v4.nb.qza \
--i-reads rep-seqs-dada.qza \
--o-classification taxonomy-dada.qza \
--p-n-jobs 5 \
--verbose
qiime metadata tabulate \
--m-input-file taxonomy-dada.qza \
--o-visualization taxonomy-dada.qzv
qiime taxa barplot \
--i-table table-dada.qza \
--i-taxonomy taxonomy-dada.qza \
--m-metadata-file metadata.tsv \
--o-visualization taxa-bar-plots-dada.qzvEn este gráfico podemos ver la variedad de especies encontradas en las muestras (Figura 7). Como vemos en la figura de la izquierda, cada muestra tiene una variedad distinta. Cuando comparamos la variedad relativa de las muestras por universidades y sexo vemos que existen poblaciones distintas según el grupo.
Figura 7: Abundancia relativa de taxones por muestra, universidad y sexo.
Para observar si hay diferencias en la composición microbiana entre ratones de las distintas universidades, aplicaremos DESeq2 a los conteos taxonómicos. Esta herramienta realiza ajustes para la profundidad de secuenciación y variabilidad biológica, y mediante modelos estadísticos, nos ayudará a identificar de manera fiable y precisa los taxones cuya abundancia difiere significativamente entre las muestras.
dds <- DESeqDataSetFromMatrix(countData=counts.taxa,
colData=metadata,
design=~University, tidy = TRUE)
dds <- DESeq(dds)
res <- results(dds)
res <- res[order(res$padj),]
rownames(res) <- sub(".*\\.(.*)$", "\\1", rownames(res))
res <- data.frame(res)Observamos que las especies más significativamente distintas entre universidades se muestran en la tabla siendo la especie Prevotella la más significativa con un log2FoldChange de 10 veces más en la Universidad de Leuven que en la Universidad de Murcia (Tabla 4).
| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| Prevotella | 115.710659 | 10.035202 | 1.075937 | 9.326945 | 1.089666e-20 | 1.155046e-18 |
| Parasutterella | 58.977861 | -25.773993 | 2.961139 | -8.704080 | 3.201618e-18 | 1.696857e-16 |
| distasonis | 60.277236 | 9.090762 | 1.081676 | 8.404334 | 4.302956e-17 | 1.520378e-15 |
| massiliensis | 5.395228 | 21.154224 | 2.972465 | 7.116728 | 1.105195e-12 | 2.928766e-11 |
| caecigallinarius | 129.895009 | 10.186081 | 1.581145 | 6.442218 | 1.177398e-10 | 2.496085e-09 |
| veroni666751 | 1459.347035 | -9.119939 | 1.439587 | -6.335106 | 2.371776e-10 | 4.190137e-09 |
En el volcano plot podemos ver las especies más representadas en una Universidad o en otra (Figura 8). Esas mismas especies podrían ser las responsables que los ratones de la Universidad de Leuven se vean menos afectados al tratamiento por LPS y no entren en un estado de sepsis severo.
Figura 8: Volcanoplot que muestra la relación entre el cambio (log2 Fold Change) y la significancia estadística (-log10 p-valor) de especies entre universidades. Las especies más presentes en la Universidad de Leuven están representados en amarillo, los de la Universidad de Murcia en morado y los no significativos en gris.
Utilizaremos un modelo murino de sepsis para investigar los mecanismos celulares y moleculares que subyacen a la neurodegeneración y la neuroinflamación. Nos centraremos en componentes del inflamasoma y emplearemos RNAseq para ampliar nuestra comprensión de la respuesta del SNC a la inflamación sistémica y sus diferencias entre géneros, según se describe en estudios previos (Rodríguez-Ramírez et al., 2024).
Inicialmente, se recibieron doce archivos en formato FASTQ, correspondientes a doce muestras distintas: seis de machos (M) y seis de hembras (F), subdivididas en controles y tratados con LPS.
Para asegurar la calidad de las secuencias, se utilizó la herramienta FASTQC, disponible en el servidor Galaxy del máster de bioinformática. FASTQC evalúa la calidad de las secuencias de datos de secuenciación de alto rendimiento, proporcionando información sobre la calidad del base calling, la distribución de la calidad de las secuencias, el contenido de bases y la presencia de secuencias contaminantes, entre otros (Figura 9A y B).
Aunque la mayoría de las muestras mostraron una calidad aceptable, las 10 primeras y últimas bases presentaron problemas de calidad en varias muestras. Por lo tanto, se decidió recortar estas regiones utilizando Fastq Trimmer, una herramienta que permite recortar secuencias de FASTQ en posiciones específicas. Posteriormente, se realizó una segunda evaluación de calidad con FASTQC para confirmar la mejora en la calidad de las secuencias (Figura 10 A y B).
Para el alineamiento de las secuencias, se utilizó HISAT2, una herramienta de alineamiento que permite el mapeo eficiente de secuencias de ARN contra un genoma de referencia. En este caso, se empleó el genoma de ratón (versión mm9.fa.gz), disponible en el servidor Galaxy. HISAT2 es particularmente útil para el alineamiento de RNAseq por su capacidad de identificar empalmes de exones y adaptarse a la variabilidad genética (Figura suplementaria 1).
Posteriormente, se evaluó la calidad del alineamiento mediante Samtools stats, una herramienta que proporciona estadísticas detalladas del alineamiento, incluyendo la proporción de lecturas correctamente alineadas, la cobertura del genoma, y la distribución de las lecturas respecto a regiones genómicas específicas (Figura suplementaria 2).
Para la cuantificación de la expresión génica, se empleó StringTie, que utiliza los archivos BAM generados por HISAT2 para ensamblar y cuantificar los transcritos. StringTie es capaz de reconstruir con precisión los transcritos y estima las abundancias de estos utilizando un modelo estadístico que permite la comparación entre muestras. Para este proceso, se utilizó el archivo de anotaciones gencode.vM10.annotation.gtf como referencia (Tabla 5).
Una vez obtenidos los conteos de expresión de cada muestra mediante StringTie, se descargaron todos los archivos de conteo. Utilicé una función personalizada en Python para concatenar estos archivos en una sola tabla de datos (ver fichero adjunto Merge_Count_tables.py), que luego se utilizó para realizar el análisis estadístico y bioinformático en R.
El análisis en R comenzó con la importación de los datos de expresión génica combinados usando la función read.csv para cargar los datos desde un archivo CSV. Se realizó una limpieza inicial de datos para asegurar que no hubieran duplicados, utilizando la función aggregate para sumar los conteos en caso de encontrar genes duplicados. Esta preparación de datos aseguraba que cada gen estuviera representado de forma única en el análisis subsiguiente.
Para una primera comprensión de los datos, se emplearon técnicas de análisis exploratorio como el Análisis de Componentes Principales (PCA), que facilita la visualización de las diferencias globales entre las muestras. Utilizando las funciones DESeq para la normalización y plotPCA de la librería DESeq2, se analizaron las principales variaciones en los datos, permitiendo identificar agrupaciones naturales o outliers.
El núcleo del análisis consistió en el uso de DESeq2, un paquete de Bioconductor para análisis de expresión diferencial que utiliza modelos estadísticos basados en la distribución negativa binomial para testar diferencias en la expresión génica. Este análisis se realizó comparando las muestras de tratados versus controles, así como entre machos y hembras, permitiendo identificar genes que mostraban cambios significativos de expresión bajo diferentes condiciones.
Para la visualización de los resultados del análisis diferencial, se utilizaron varios enfoques:
Para interpretar los resultados en un contexto biológico más amplio, se realizaron análisis de enriquecimiento utilizando clusterProfiler y org.Mm.eg.db para asociar los genes diferencialmente expresados con rutas biológicas conocidas. Esto incluyó el uso de la función enrichGO para el enriquecimiento de términos de ontología génica y enrichKEGG para identificar rutas metabólicas clave implicadas en las respuestas observadas. Los resultados de estos análisis se visualizaron utilizando enrichplot y pathview para proporcionar interpretaciones gráficas de las rutas y procesos biológicos afectados.
La comparación inicial se realizó entre los machos controles y los machos tratados con LPS, con el objetivo de identificar diferencias en la expresión génica que podrían elucidar las respuestas del SNC a la inflamación sistémica.
En el análisis de PCA, observamos una anomalía con la muestra de un ratón tratado (2MTreated), que mostró una proximidad notable a las muestras de control (Figura 11A). Esto sugiere una posible falta de respuesta al tratamiento o un error técnico durante la administración del mismo. Al examinar la expresión de genes que mostraron una diferenciación clara en otros ratones tratados, confirmamos que la muestra 2MTreated no seguía el mismo patrón, destacando en particular en el gen Meg1 (maternal expressed gene 3) (Figura 11B). Por lo tanto, se decidió excluir esta muestra del análisis posterior (Figura 11C).
Una vez ajustado el conjunto de datos, realizamos representaciones gráficas como el volcano plot y el heatmap (Figura 12A y 12B). Estas visualizaciones confirmaron la presencia de múltiples genes diferencialmente expresados entre los machos tratados y los controles, con un patrón de expresión homogéneo en las muestras restantes. Los genes más destacados incluyen Cebpd (CCAAT/enhancer binding protein delta), Lrriq4 (leucine-rich repeats and IQ motif containing 4), Tmem256 (transmembrane protein 256), Ghdc (GH3 domain containing), Meg3 y un pseudogén relacionado con Meg3, todos con cambios significativos de expresión (Figura 12C).
En cuanto a las vías biológicas afectadas, el análisis de Gene Ontology señaló la producción de citoquinas por leucocitos mieloides, la producción y regulación de interleucina 1, y la regulación negativa del crecimiento celular o la apoptosis como procesos significativamente alterados. Adicionalmente, el análisis de vías de KEGG identificó la Enfermedad Hepática Alcohólica, la vía de señalización de ErbB, la formación de trampas extracelulares de neutrófilos y la señalización de HIF-1 como rutas destacadas (Figura 13A, B y C).
Estos resultados sugieren que la respuesta a LPS en ratones machos puede estar mediada por alteraciones significativas en la expresión de genes clave y la activación de vías moleculares críticas asociadas con la respuesta inflamatoria y la regulación celular, proporcionando así una base para futuras investigaciones sobre los mecanismos subyacentes de la neuroinflamación y neurodegeneración en contextos de sepsis.
En nuestro estudio sobre la respuesta inflamatoria inducida por LPS en ratones, la comparación entre las hembras control y las hembras tratadas con LPS reveló hallazgos significativos, destacando diferencias claras en la expresión génica. Utilizando el análisis de PCA, confirmamos la validez de todas las muestras incluidas en el estudio, sin la necesidad de excluir ninguna por comportamiento atípico (Figura 14A).
El análisis mediante volcano plot demostró un número superior de genes diferencialmente expresados en las hembras en comparación con los machos, acompañado de mayores valores de Fold Change y p-values más bajos, indicando cambios más pronunciados y estadísticamente más significativos (Figura 15A). El heatmap reflejó un comportamiento consistente entre las muestras de ratones control y tratados, con patrones de expresión uniformes para los genes representados (Figura 15B).
Dentro de los seis genes más diferencialmente expresados, encontramos coincidencias con los machos en genes como Meg3, Lrriq4 y Ghdc. Sin embargo, en las hembras también se identificaron genes no anotados como Gm43473 y Gm35339, junto con Mcm10 (minichromosome maintenance 10 replication initiation factor), que son de particular interés debido a sus potenciales roles en la replicación del ADN y la regulación celular (Figura 16).
Respecto a las vías biológicas alteradas, observamos la misma actividad en las rutas identificadas en los machos, con la adición significativa de la ruta de regulación de la contracción del músculo liso en las hembras (Figura 17A y B). El análisis de las vías de KEGG resaltó la participación de vías implicadas en varios tipos de cáncer, como el cáncer de páncreas y la leucemia mieloide, así como procesos metabólicos específicos como el metabolismo de la colina y la ruta de señalización AGE-RAGE (Figura 17C).
Es crucial destacar que, en las hembras, se observó la desregulación de 254 genes más que en los machos, lo que refleja una mayor sensibilidad a la sepsis inducida por LPS (Figura 18). Este fenómeno subraya la importancia de considerar las diferencias de sexo en la respuesta biológica a tratamientos y condiciones patológicas. La identificación detallada de estos genes proporciona una base sólida para investigar los mecanismos subyacentes a estas diferencias de respuesta y podría facilitar el desarrollo de estrategias terapéuticas más personalizadas y efectivas, tomando en cuenta las particularidades biológicas de cada sexo biológico.
Mediante el análisis de volcano plot y heatmap, hemos identificado una serie de genes que exhiben una expresión diferencial significativa entre machos y hembras (Figura 19A y B).
Entre los genes más destacados se encuentran Lima1 (LIM Domain And Actin Binding 1), que está implicado en la reorganización del citoesqueleto; Kif4 (Kinesin Family Member 4A), que juega un papel en el transporte intracelular y la segregación de cromosomas durante la división celular; Trpc4ap (Transient Receptor Potential Cation Channel Subfamily C Member 4 Associated Protein), relacionado con la señalización de calcio; Kpna1 (Karyopherin Subunit Alpha 1), que está involucrado en el transporte nuclear; Ano7 (Anoctamin 7), implicado en procesos de señalización celular y Cdc6 (Cell Division Cycle 6), esencial para la iniciación de la replicación del ADN (Figura 20).
Un análisis complementario mediante un diagrama de Venn permite cuantificar la distribución de los genes diferencialmente expresados entre machos y hembras enfermos. Se identificaron 1134 genes exclusivos en machos, lo que representa el 23.3% del total de genes diferencialmente expresados, mientras que en las hembras se encontraron 1391 genes únicos, constituyendo el 28.5% del total. La intersección entre ambos sexos comprende 2351 genes, representando el 48.2% del total, evidenciando una respuesta común considerable al tratamiento (Figura 21). Este dato resalta la existencia de una base común en la respuesta al LPS, aunque también se observan diferencias significativas que reflejan respuestas específicas de cada sexo.
Este patrón de expresión refleja que, si bien la principal influencia en la expresión génica parece ser el tratamiento con LPS, las diferencias entre machos y hembras también son notables, sugiriendo que la variabilidad genética asociada al sexo contribuye a modulaciones específicas en la respuesta molecular. Esto es particularmente relevante en el contexto de una medicina personalizada, donde comprender estas diferencias puede permitir la optimización de tratamientos basados en características genéticas individuales, incluyendo el sexo del paciente.
En nuestro estudio, el análisis scRNAseq se centró exclusivamente en modelos murinos machos, a fin de desentrañar las respuestas celulares específicas involucradas en la sepsis. Para ello, se aislaron células positivas para CD11b (indicativas de macrófagos y microglia) y CD45 (marcador de células linfoides), utilizando estos marcadores para diferenciar el estado de activación y el tipo celular dentro del contexto inflamatorio de la sepsis. En este caso, las muestras a comparar corresponden a microglía extraída de ratones control y ratones tratados con LPS, todos machos.
La selección de células positivas tanto para CD11b como para CD45 nos permitió enfocarnos en macrófagos de microglia activados, células que desempeñan roles cruciales en la respuesta inmunitaria y en la mediación de la neuroinflamación dentro del SNC. Estos macrófagos activados son indicativos de una respuesta inmune intensificada y juegan un papel fundamental en la patogénesis de la sepsis, facilitando tanto la defensa contra patógenos como la potencial contribución al daño tisular asociado.
Por otro lado, las células que presentaban positividad exclusiva para CD11b fueron consideradas como microglia en un estado de activación menos pronunciado. Estas células, aunque participan en la respuesta inmune, se encuentran en un estado de activación menos agresivo comparado con los macrófagos microgliales activados, lo que sugiere una respuesta inmunitaria más modulada o en etapas iniciales de activación.
Finalmente, las células exclusivamente positivas para CD45 se interpretaron como componentes restantes del compartimento linfoides, englobando una variedad de células inmunitarias implicadas en la respuesta a la sepsis. Esta fracción representa una gama diversa de células inmunológicas, incluyendo linfocitos T, B, y otras células linfoides, que participan en la coordinación y ejecución de la respuesta inmunitaria adaptativa y que pueden jugar un papel en la progresión de la neuroinflamación y neurodegeneración en el contexto de la sepsis.
Antes de iniciar el análisis, se requiere la configuración adecuada del entorno de trabajo y la preparación de los datos iniciales. El análisis comienza con archivos FASTQ procesados por la plataforma de bioinformática del IMIB utilizando el software Cell Ranger, un método estándar en la industria para el procesamiento de datos de secuenciación de células únicas obtenidos mediante tecnologías como 10x Genomics. Los archivos resultantes, específicamente barcodes.tsv, features.tsv y matrix.mtx.gz, contienen la información esencial para la creación de objetos de análisis y son el punto de partida de nuestro estudio.
La primera etapa del análisis en R incluye la configuración del entorno y la carga de paquetes necesarios. Se utilizan diversas bibliotecas de R, que facilitan desde la manipulación de datos hasta visualizaciones avanzadas. Entre las más destacadas están: - Seurat: esencial para análisis de secuenciación de células únicas. - dplyr, ggplot2, tidyverse: para manipulación y visualización de datos. - patchwork, cowplot: para la combinación y organización de gráficos. - sctransform, ReactomeGSA, topGO: para normalización y análisis de enriquecimiento. - SingleCellExperiment, SingleR, celldex: para manejo de experimentos de células únicas y anotación de tipos celulares.
Utilizando la función Read10X de Seurat, se cargan los datos desde los directorios especificados. Para cada conjunto de datos (Control y LPS), se crean objetos Seurat a partir de las matrices de expresión. Estos objetos son fundamentales para el manejo eficiente de los datos y facilitan la realización de análisis subsiguientes.
El análisis comienza con la evaluación de la calidad de las células, utilizando métricas como el número de genes detectados y el porcentaje de ADN mitocondrial. Se aplican filtros para eliminar células de baja calidad que puedan afectar los resultados del análisis.
Los datos son normalizados para eliminar variaciones técnicas. Posteriormente, se identifican características (genes) que muestran alta variabilidad, lo cual es crucial para las etapas de reducción de dimensionalidad y agrupación.
Se utilizan técnicas como PCA (Análisis de Componentes Principales) para reducir la dimensionalidad de los datos. Este paso es seguido por métodos de agrupación para identificar grupos o clústeres de células con perfiles de expresión similares.
El análisis culmina con la integración de todos los datos y resultados. Se realizan visualizaciones avanzadas como UMAP y t-SNE para representar las relaciones entre las células y clústeres. Además, se utilizan gráficos como heatmaps y dotplots para representar los datos de una manera que facilite la interpretación y presentación de los resultados.
Para cada clúster identificado, se buscan genes que se expresan de manera única, lo que ayuda a caracterizar los tipos celulares dentro de los clústeres.
Se comparan las células bajo diferentes condiciones experimentales (Control vs. LPS) para identificar cambios significativos en la expresión génica que puedan ser indicativos de respuestas biológicas al tratamiento.
Para los genes de interés, se realizan análisis de enriquecimiento utilizando herramientas como ReactomeGSA y topGO para identificar rutas biológicas o procesos celulares relevantes.
Los resultados obtenidos en el análisis de single-cell RNA-seq revelan aspectos cruciales sobre el impacto biológico del tratamiento con LPS en ratones, principalmente evidenciando una activación significativa del sistema inmunitario y cambios en la composición celular que podrían ser indicativos de una respuesta inmunitaria compleja o incluso de procesos patológicos similares a los observados en enfermedades neurodegenerativas.
Una de las principales limitaciones observadas fue el bajo número de células aisladas para el experimento, lo que podría afectar la representatividad de las poblaciones celulares y la fiabilidad de los resultados en términos de expresión diferencial de genes (Figura 22). Esta limitación es crucial, pues podría conducir a una sobrestimación de ciertos efectos o a la no detección de poblaciones celulares relevantes.
El análisis permitió identificar siete clústeres celulares principales, incluyendo oligodendrocitos, astrocitos, células endoteliales, monocitos, macrófagos activados, neuronas, y otras células del sistema inmunitario como células NK, células B y células T (Figura 23A). La presencia aumentada de monocitos y células del sistema inmunitario en el ratón tratado con LPS sugiere una movilización de estas células en respuesta a la sepsis. Asimismo, el incremento de oligodendrocitos podría estar vinculado a la sobreexpresión de rutas relacionadas con la diferenciación de células gliales observada en análisis anteriores (Figura 23B).
La sobreexpresión de genes en células inmunitarias y microglía, particularmente en condiciones de estimulación con LPS, destaca genes como Atp5c1 o H2afz (Figura 24). Curiosamente, todos los genes diferencialmente expresados mostraron sobreexpresión, lo cual refuerza la idea de una activación aguda del sistema inmunitario y de respuestas celulares en el contexto de inflamación y posible daño tisular.
Las rutas identificadas a través de bases de datos como Reactome, Gene Ontology (GO) y KEGG reflejan una amplia gama de procesos biológicos afectados:
Catozzi C, Cuscó A, Lecchi C, Carlo ED, Vecchio D, Martucciello A, D’Angelo L, Francino O, Bonastre AS, Ceciliani F (2019) Impact of intramammary inoculation of inactivated Lactobacillus rhamnosus and antibiotics on the milk microbiota of water buffalo with subclinical mastitis. PLoS ONE 14(1): e0210204. <https://doi: 10.1371/journal.pone.0210204>
Rodríguez-Ramírez, K. T., Norte-Muñoz, M., Lucas-Ruiz, F., Gallego-Ortega, A., Calzaferri, F., García-Bernal, D., Martínez, C. M., Galindo-Romero, C., de Los Ríos, C., Vidal-Sanz, M., & Agudo-Barriuso, M. (2024). Retinal response to systemic inflammation differs between sexes and neurons. Frontiers in immunology, 15, 1340013. <https://doi:10.3389/fimmu.2024.1340013>
Salipante SJ, Kawashima T, Rosenthal C, Hoogestraat DR, Cummings LA, Sengupta DJ, Harkins TT, Cookson BT, Hoffman NG. Performance comparison of Illumina and ion torrent next-generation sequencing platforms for 16S rRNA-based bacterial community profiling. Appl Environ Microbiol. 2014 Dec;80(24):7583-91. doi: 10.1128/AEM.02206-14. Epub 2014 Sep 26. Erratum in: Appl Environ Microbiol. 2016 Jul 15;82(14):4453. PMID: 25261520; PMCID: PMC4249215.